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When a blood vessel ruptures or gets inflamed, the 
human body responds by rapidly forming a clot 
to restrict the loss of blood. Platelets aggregation 
at the injury site of the blood vessel occurring via 
platelet-platelet adhesion, tethering and rolling on the 
injured endothelium is a critical initial step in blood 
clot formation. A novel three-dimensional multi-scale 
model is introduced and used in this paper to simulate 
receptor-mediated adhesion of deformable platelets 
at the site of vascular injury under different shear 
rates of blood flow. The novelty of the model is based 
on a new approach of coupling submodels at three 
biological scales crucial for the early clot formation: 
novel hybrid cell membrane submodel to represent 
physiological elastic properties of a platelet, stochastic 
receptor-ligand binding submodel to describe cell 
adhesion kinetics and lattice Boltzmann submodel for 
simulating blood flow. The model implementation on 
the GPU cluster significantly improved simulation 
performance. Predictive model simulations revealed 
that platelet deformation, interactions between 
platelets in the vicinity of the vessel wall as well as 
the number of functional GPIba platelet receptors 
played significant roles in platelet adhesion to the 
injury site. Variation of the number of functional 
GPIba platelet receptors as well as changes of 
platelet stiffness can represent effects of specific drugs 
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reducing or enhancing platelet activity. Therefore, predictive simulations can improve the 
search for new drug targets and help to make treatment of thrombosis patient-specific. 



1. Introduction 

When a blood vessel ruptures or gets inflamed, the human body responds by rapidly forming a 
clot to restrict the loss of blood. Blood clot (thrombi) formation is a complex biological process 
involving an extensive system of biochemical (coagulation) reactions, platelet hydrodynamics, 
platelet-platelet and platelet-blood vessel wall interactions leading to ligand-receptor adhesion 
bond formation and platelet activation. Platelet adhesion to the vessel wall is one of the first 
events associated with formation of haemostatic clots and pathological thrombi. 

In this paper, a new three-dimensional multi-scale model of platelet-blood flow-vessel 
wall interactions combining submodels at three biological scales crucial for the early platelet 
aggregation is introduced and calibrated to investigate how platelet stiffness, GPIb receptor 
expression and platelet-platelet interaction affect platelet-wall adhesion quantified in terms of 
platelet pause time. We implemented a novel approach of combining a recently developed platelet 
hybrid membrane submodel, the subcellular element (SCE) representation of the cytoskeleton 
network and a continuum description of the lipid bilayer to study the very first step of blood clot 
formation, the rapid formation of unstable bonds which slow platelets and cause platelet-flipping 
and adhesion to the damaged surface. The hybrid platelet model was also coupled with the 
lattice Boltzmann model (LBM) of blood using the immersed boundary method (IBM) to simulate 
platelet motion and deformation in shear flow. The kinetics-based adhesive dynamics model 
was also integrated into the three-dimensional model to simulate formation and disassociation 
of the receptor-ligand bonds during the platelet-platelet and platelet-vessel wall interactions. 
Parallelized model simulations were implemented on a GPU computer cluster which speeded up 
simulations by a factor of 100 (table 3) in comparison with CPU implementation which allowed, 
for the first time, to run biologically relevant predictive simulations. 

By using a novel biologically calibrated three-dimensional modelling approach, it is shown 
that the platelet stiffness, the number of GPIba platelet functional receptors and mutual 
interaction between platelets can significantly alter the adherence of platelets at the site of 
vascular injury. Our results demonstrate how a comprehensive modelling approach coupling 
three biologically relevant scales can provide new insights into the biomedically important 
problem of early thrombus development. Variation of the number of functional GPIba platelet 
receptors as well as changes of platelet stiffness can represent effects of specific drugs for reducing 
or enhancing platelet activity. This emphasizes the importance of predictive simulations as they 
can potentially improve the search for new drug targets and help with making treatment of 
thrombosis patient-specific. 

Damage or alteration of a blood vessel lining can result in activation of flowing platelets 
and their subsequent aggregation at sites of vascular injury. The ability of platelets to tether to 
and translocate on injured vascular endothelium relies on the interaction between the platelet 
glycoprotein receptor Iba(GPIba) and the Al domain of von Willebrand factor (vWF-Al) [1]. 

Along with biochemical activation of platelets, large shear disturbances of blood flow are 
one of the key factors promoting pathogenic activation of platelets and formation of thrombi. 
In addition, platelets flowing in a whole blood exhibit increased concentrations in the vicinity 
of the vessel wall, making platelet-platelet interactions more frequent near vascular surfaces. 
Excessive accumulation of platelets at injury sites is one of the pathological events that result 
in acute myocardial infarction, sudden death and ischaemic stroke. This pathological process is 
responsible for mortality and morbidity rates higher than for any other disease, making platelets 
a major target for therapeutic interventions. Thus, studying individual platelet dynamics as well 
as platelet-platelet interactions and platelet adhesion to a vascular or thrombus surface is of high 
biomedical importance and urgency. 



Effects of shear flow on accumulation of platelets on various surfaces have been extensively 
studied in in vitro and in vivo experiments [1-7]. However, there is a limited amount of available 
experimental data on an individual platelet dynamics in the vicinity of the vascular surface as 
well as platelet-surface attachment. There is also a lack of experimental data demonstrating 
how platelet-surface attachment is affected by mechanical properties of a platelet as well as by 
platelet receptor-ligand kinetics. Better understanding of platelet aggregation requires study of 
the interplay among biochemical, mechanical and hydrodynamic processes occurring at different 
scales, including a nanometre scale (receptor-ligand kinetics), a micrometre scale (cellular level) 
and a millimetre scale (early platelet aggregate). Multiple characteristic scales make it difficult 
to experimentally discern effects of different processes involved in platelet-surface attachment 
and overall thrombus growth dynamics. Meanwhile, a multi-scale modelling approach can 
provide a useful predictive tool to aid in elucidating mechanisms of platelet-wall attachment 
and aggregation. 

Several multi-scale models attempting to couple large numbers of submodels at different 
scales have been developed (see, among others, for reviews [8,9]). These models implemented 
simplified submodels in order to make simulations less computationally expensive. It is extremely 
difficult at this time, if not impossible, to validate predictions of multi-scale models attempting 
to combine submodels at all scales representing processes of blood clot formation using existing 
experimental data. In addition, most experimental data are currently available at the molecular 
level and individual platelet level. Therefore, it is important to develop detailed multi-scale 
models coupling two or three scales and considering only a few processes at a time. Such 
models when properly calibrated with available experimental data can provide useful predictive 
tools aiding in designing new experiments, drug design and planning new patient-specific 
therapeutic strategies. 

Several computational models have been developed to characterize platelet and other types of 
blood cell motion and adhesion dynamics under hydrodynamic shear flow at cell and receptor 
levels (see [8,9] for a review). Analytical solutions for forces and torques exerted on a platelet 
treated as a rigid object in the Stokes regime in a two-dimensional case were obtained in 
[1,10] and compared with the data obtained using an image analysis algorithm for tracking 
the motion of platelets before, during and after contact with the surface. Kinetic properties of 
the receptor-ligand adhesion bonds, GPIba-vWF, were quantified in [1,4] using Monte Carlo 
simulations and pause time analysis of transient capture/ release events. This approach provided 
association and disassociation rate constants k on and k Q ft, depending on the shear rate of the 
blood flow. 

Experimental study in [11] showed that platelets have viscoelastic properties and the elastic 
moduli in the range of 1-50 kPa. Large deformation occurred when platelets were suspended in 
the shear flow [12]. To account for the elastic and viscoelastic properties of cells, a number of 
methods accounting for cell structural properties have been developed [13,14]. The SCE model 
introduced in Sandersius & Newman [13] represented each cell by a collection of elastically 
coupled SCEs, interacting with each other via short-range potentials. Sweet et al. [15] and Xu 
et al. [14] presented a three-dimensional modelling approach in which cells, modelled by SCEs, 
were coupled with fluid flow and substrate models by using the Langevin equation. 

The fluid-structure interaction approach is an essential part of the model. Previously, the IBM 
introduced by Peskin [16] to investigate blood flow in the human heart has been applied to many 
other fluid-structure interaction problems, including platelet aggregation [17] and deformation 
of red blood cells [18]. Skorczewski et al. [19] developed a two-dimensional model, using a 
lattice Boltzmann IBM to investigate the motion of platelets near a vessel wall and close to an 
intravascular thrombus, in which they modelled the platelets as rigid bodies, whereas the red 
blood cells were represented as deformable bodies. 

The results of the predictive simulations of the three-dimensional model introduced in this 
paper revealed that the platelet pause time strongly depends on the stiffness of the platelet 
as well as on the number of expressed GPIb membrane functional receptors. Additionally, we 
demonstrated that the platelet-platelet interaction near the surface of the vessel wall could 



significantly decrease the platelet paused time, and thus decrease the rate of platelet attachment 
to the injury site. 

This paper is organized as follows. It starts with the Biological Background section. Then the 
Methodological Innovation is described in detail, including description of submodel at each of 
three space scales and of the coupling approach. This is followed by the Results section, which 
includes model validation and description of the predictive simulation. Biological relevance of 
the predictive simulations is discussed in the Discussion section. GPU implementation of the 
three-dimensional model is described in appendix A. 

2. Biological background 

The mechanism by which platelets bind to a damaged blood vessel wall is similar to that of 
leucocyte binding to activated endothelium [20], and requires two binding steps. The first step 
is the rapid formation of unstable catch-slip bonds which slow platelets and cause platelet- 
flipping along the damaged surface. (Counterintuitively, the dissociation rate first decreases 
with increasing force until reaching a threshold.) This is mediated by the platelet receptor 
component, GPIba, forming transient bonds with the vWF exposed at the injury site. Rapid 
association and dissociation kinetics of the bonds result in transient tethering and subsequent 
flipping (or rolling) and pausing of platelets on the vessel surface [1,21]. Then, stable bonds 
slowly form between platelet receptors and ligands (often integrin aTIb/33 binding with vWF 
or fibrinogen) bound to the damaged wall or the surface of the thrombus resulting in strong 
adhesion, initiating transmembrane and, subsequently, intracellular signalling. As the blood clot 
grows, platelet-platelet interaction becomes one of the major factors determining clot growth rate 
and integrity as platelets expose GPIIbllla receptors which permit platelet-platelet adhesion via 
fibrinogen. Adhesion of platelets to the injured surface is also affected by shear rates of the flow. 
At high shear, platelet integrin a2/31 and GPVI receptors are not sufficient to initiate binding to 
collagen, and binding of the GPIba receptor to vWF immobilized on collagen becomes essential in 
platelet adhesion. 

The stiffness of the platelet not only determines the shape and morphology of the clot, 
but also affects clot mechanical properties, as platelet stiffness determines cell shape when 
exposed to various flow conditions and contact interaction with other cells and blood vessel wall. 
This will affect the number of receptor-ligand pairs in platelet-platelet and platelet-substrate 
interactions. Platelet stiffness is also an important property reflecting platelet functioning, 
because it reorganizes its structure during activation or as a response to physiological or 
pathological conditions. 

To date, to the best of our knowledge, cumulative effects of platelet stiffness, different levels of 
expression of GPIba receptors and platelet-platelet interaction impacting the strength of platelet- 
substrate binding have not been systematically investigated. Our model provides a unique means 
for quantitatively understanding these effects, which are critical for improving our knowledge 
about the initial stage of the blood clot formation. 

3. Methodological innovation of the three-dimensional modelling approach 

The novelty of the three-dimensional model lies in developing a novel membrane submodel as 
well as in new approaches of coupling submodels of biological processes at three spatial scales 
(figure 1) which are crucial to early blood clot formation. At the subcellular scale (nanoscale), a 
kinetics-based stochastic dynamic adhesion submodel (DAM) is used to simulate vWF-GPIba 
binding and GPIba-vWF-GPIba binding, in which individual vWF and GPIba molecules are 
represented as elastic springs (table 1). This is justified by the fact that these receptor-ligand 
bindings are probabilistic in nature [1]. Moreover, individual filaments in the cytoskeleton 
network of the platelet membrane are treated as coarse-grained harmonic springs. At the cellular 
scale, a novel continuum description of the lipid bilayer of the cell membrane is used. We 





Figure 1. Flow chart of the simulation algorithm. Using Fj(t) hybrid membrane model, the forces Fj(t) acting on cell elements 
were calculated. The forces acting on fluid nodes were spread from Fj(t) by the immersed boundary method. The velocity 
field v[x,t + Af) of fluid was obtained by the lattice Boltzmann method. The velocities of cell elements d>Y/df(f + Af)were 
determined by immersed boundary interpolation. The cell elements positions were updated based on the velocities. Finally, the 
stochastic adhesion model was used to determine the force fbond acting on receptor-ligand bonds that bind platelets to vessel 
walls. (Online version in colour.) 



Table 1. Biological processes and submodels at different scales. 



<0.1 ixmsubcellular-level 
nanoscale 



processes 



ligand-receptor 
interactions 



approximately 1 |xm 
cellular-level microscale 



individual platelet 
deforming, flipping 
and adhering to 
vessel wall 



>10 |xm macroscale 



blood flow and its 
Interaction with 
platelet 



submodels 



stochastic dynamic 
adhesion model 



hybrid membrane 
model 



lattice Boltzmann 
method 



coupling 



1. Nano-micro scales: coupled by 

explicitly modelling receptors 
on platelet membrane nodes 

2. Micro-macro scales: coupled 
through immersed boundary 
method 



developed this new platelet membrane model to study effects of membrane stiffness on cell- 
substrate interaction, which was shown to strongly affect platelet-injury site adhesion. (See also 
§4c(i) for model prediction.) The subcellular scale and the cellular scale components are integrated 
by distributing GPIba receptors at the vertices of the cytoskeleton network and by superimposing 
the cytoskeleton network and the lipid bilayer. At the macroscale, the dynamics of the fluid 
flow is represented using the LBM to facilitate parallelizing the simulation code on GPUs. The 
platelet model is coupled with the LBM using the IBM. (The coupling and data flow between all 
the submodels are demonstrated in figure 1.) We calibrate and validate this three-dimensional 
model by comparing simulations at different scales with either theoretical results or available 
experimental data at these scales. Specifically, the platelet model coupling with the LBM was 
validated using theoretical results and previous simulation results (see also §4a), whereas the 
platelet-substrate adhesion simulations were compared with experimental data to calibrate the 
DAM under different flow conditions (see also §4fr). 

At each time step of simulation, the hybrid membrane model is first used to calculate forces 
acting on the nodes of the Lagrangian mesh representing platelet geometry, such as bond 
forces resulting from stretching or compression of the cytoskeleton network, bending forces 
resulting from deformation of the lipid bilayer and attraction/ repulsion between platelets and the 
environment owing to formed ligand-receptor bonds. This is followed by coupled LBM and IBM 
to update fluid flow and position of platelets. Finally, the MC computations of platelet adhesion 
to a surface expressing vWFs are performed to break the already formed bonds and to generate 
new bonds from unbound GPIba and vWF. 

We note that this is the first time that a detailed platelet membrane model has been developed 
and implemented on GPUs for studying cell-flow, cell-cell and cell-substrate interactions. 
Because of the speed-up gained by GPU implementation, we are able to investigate effects of 
these interactions and cell mechanics on platelet dynamics in a timely manner. Additionally, this 
model can be directly used for modelling any biological cells with membrane structures similar 
to these of eukaryotic cells. Here, we describe in detail individual submodels and explain how 
they are coupled. 

(a) Platelet membrane submodel 

We simulate the motion of platelets in a three-dimensional region bounded by an infinite flat 
plane at z = 0 (see figure 2a for example). A platelet has initial shape defined by x 2 /a 2 + y 2 /a 2 + 
z 2 /(Xa) 2 = 1, where a = 1 |xm is the approximate particle radius and X = 0.25 is the aspect ratio 
[22,23]. The Reynolds number of this system is Re = ypa 2 /fi = 0(1(T 3 ), where y = 300 and 400 s" 1 
are the shear rates used in experiments [1], a = 1 |xm is the particle radius, p = 1.0239 g cm -3 is the 
density of blood plasma and \i — 1.2 cP is the viscosity of blood plasma [24]. 

The platelet membrane, which is similar to the membrane of a red blood cell, is also assumed to 
consists of a lipid bilayer and an attached cytoskeleton. Following ideas from [16,25], the platelet 
membrane surface geometry is represented by a triangular mesh consisting of a collection of N 
(N = 958 in our simulation) points {X z -, iel...N) (figure 2b). The connected edges of the mesh are 
used to model the cytoskeleton network of the platelet membrane, and the triangulated mesh 
surface represents the lipid bilayer of the cell membrane, where the cytoskeleton attaches to. 
The mesh points represent coarse-grained actin vertices and each edge of the mesh represents a 
coarse-grained filament. The Helmholtz free energy of the membrane is defined as 

^membrane — ^SCE + ^bending + ^volume + ^area + ^wall- (3-1) 

Here, Hsce is the in-plane energy of the cytoskeleton network; H^ en ^ ins is the bending energy 
representing the bending resistance of the lipid bilayer; H vo i ume and H area are volume and area 
conservation constraints, respectively, and H wa u represents the energy relating to interactions due 
to ligand-receptor binding (explained in detail in §3c). 

We use a harmonic 'spring' model to simulate the elasticity of the edge connecting mesh points 
i and j, which mimics a coarse-grained filament. The associated potential energy functions for 




Figure 2. (a) Schematic diagram showing one platelet translating and rotating in shear flow near an infinite plane wall. 
(b) Structure of a platelet consisting of 958 SCEs. The major radius, a, and centroid height, H, are defined as is the coordinate 
system and flow direction. One platelet is represented by a collection of elastically linked SCEs, interacting with one another via 
a spring-like elastic force. The GPIba receptors are randomly uniform distributed on the platelet membrane, and vWF ligands 
are distributed on the wall. (Online version in colour.) 



points i and j are 



k / 



LT|=-(||R iy ||-L^, (3.2) 

where Ljj is the rest length, Rjj = Xj — X\ the position vector difference for points i and /', 
respectively, and k = EAx/5 the coefficient that defines the spring 'stiffness' [26] for elastic 
modulus E = 25 kPa [11] and Ax = 0.1 |xm unit link length of the spring. The total potential energy 
for the cytoskeleton network is Hsce = for all edges The corresponding force vector acting 
on point i by point j is 

f?. = - V x U?.(x) = -fc(||R, y || - 1^ ) A. (3.3) 

The area and volume conservation constraints, which account for area incompressibility of the 
lipid bilayer and incompressibility of the inner cytosol, respectively, are expressed as 

MS^-s^y _ MS -So) 2 

= — — „ t L , 

U all triangles 

and ^ 

^volume — ZVq ' P'^) 

where k S/ kt and k v are the global area, local area and volume constraint coefficients, respectively. 
The terms S total and V denote the surface area and volume for the whole platelet, whereas So, 



S* otal and V 0 

are the individual triangle mesh area, the total membrane area and the volume for 
unstressed platelet, respectively. 

We adopt the energetic variational approach developed in [27] to represent the lipid bilayer of 
the cell membrane. Let £ e R 3 be a smooth, closed surface representing the lipid bilayer of the 
platelet. The bending energy of the lipid bilayer is defined as [27] 

1 



^bending — ^0 



K{x) 2 dS{x), (3.6) 
s z 



b 



where K(x) = l/2(/q(x) + K2(x)) is the mean curvature, and k\(x), K2(x) are the principle curvatures 
at the point x. We follow the finite-element method in [28] to calculate k\(x) and K2(x). Briefly, let 
w(£, rj) be a function defined over a triangle of the surface mesh representing the lipid bilayer and 
approximated as 

6 

u(^ = ^N^), (3.7) 

z=l 

where %, r\ are the local parametric coordinates, u\ is the value of u at node i and N z - (£,??) 
are the basis functions for a quadratic six-node triangular finite element. To evaluate the 
membrane curvature tensor k, one needs to calculate the left Cauchy-Green strain tensor, which 
is determined from the surface deformation gradient tensor, A. For each triangular element, the 
surface deformation gradient tensors at the element nodes are obtained by solving the following 
system of equations ■ S 

ax ax ax ax , nns 

A = — , A = — and A h = 0, (3.8) 

3£ 3£ drj 3rj v } 

at each node of the element, and X, X are its positions in the unstressed state and after 
deformation at time t, respectively, h is the unit normal vector to the undeformed membranes. 
To evaluate the curvature tensor k at a point, one needs to solve 

ax dn ax dn , v 

— • k = — , — • k — — and n • k = 0, (3.9) 

3§ a£ dn dr) v } 

at each element node and then average over the elements sharing that node, and n is the unit 
normal vector to the deformed membranes. The mean curvature is given by 

K(x)=^(K 1 +K 2 )=hr( l c). (3.10) 

The normal component of the elastic force-associated bending energy (3.6) is obtained by taking 
variational derivative and is given by Fbend — (AzK + 2K 3 )n. 

Thus, nodal forces Fj are derived from the total energy as follows 

a(H S CE + ^volume + #area + #wall) . r /o 1 1 \ 

Fi= h Fbend- (3-11) 

dX{ 

Computation of d(H wa n)/dXj is explained in §3c. 

(b) Platelet stochastic dynamic adhesion submodel 

The kinetics-based stochastic DAM based on ideas of the Dembo model [29,30] is used to simulate 
the GPIba of unactivated platelet binding to immobilized vWF on the vessel-wall or platelet- 
platelet adhesion through forming GPIba-vWF-GPIba bonds, in which vWF was originally in 
plasma. Here, we provide details of the model for GPIba-(immobilized) vWF bond formation; 
modelling of GPIba-vWF-GPIba is treated similarly. Each platelet has approximately 10688 
GPIba receptors distributed uniformly on its membrane surface, to achieve a surface density 
of approximately 1500 receptors |xm -2 [31]. In our model, 5344 receptor point locations on 
the platelet surface are randomly distributed on the platelet membrane mesh, with each point 
location representing two GPIba receptors, because there are two GPIba receptors existing on 
each GPV molecule [32]. On the bottom plane of the simulation domain, z = 0, immobilized 



vWFs are uniformly distributed, resulting in a vWF density of 25 |xm -2 , which is consistent with 
experimental conditions in Doggett etal. [1]. 

The following rules are used for governing the GPIba-vWF binding [33]. (i) Two vWF 
molecules cannot bind to the same receptor nodes for reasons of steric blocking, and (ii) receptors 
from a maximum of four receptor nodes present on a platelet surface can bind a vWF molecule. 

In our stochastic DAM, when an unbound GPIba and an unbound vWF are separated by less 
than the length of a GPIba-vWF bond of 128 nm [34,35], a test for forming a bond is performed. 
Next, the formed bonds are tested for breakage. A GPIba-vWF bond is modelled as a linear 
spring. 

Probabilities of GPIba-vWF bond formation and dissociation are calculated using Pf 
(probability of forward reaction) and P r (probability of reverse reaction) described in [36]: 
Pf = 1 — exp(— k on At),P r = 1 — exp(— k 0 # At), where Zc 0 ff and k on are given in s _1 units and At 
is the simulation timestep. The reverse rate constant is calculated using the Bell model for the 
force-dependent dissociation rate of weak non-covalent bonds 

fc °ff=Cexp(0), (3-12) 

where fc 0 ff(£b) is the bond dissociation rate, k® {{ is the unstressed off-rate, y is the reactive 
compliance, F\, is the applied force on the bond and k^T is the product of Boltzmann constant 
and temperature. The dependence of bond formation rate constant k on on the deviation bond 
length is described by [29,33] as 

fc on = Cexp (a\x h - ; b | K-0-5|x b -i b | ^ ^ 

where k® n is the intrinsic cross-linking formation rate constant, a is the spring constant, l\> is the 
equilibrium bond length, x^ is the distance spanning the endpoints of the GPIba receptor on the 
platelet surface and the vWF-Al binding site. 

The adhesion force of the GPIba-vWF bond located at zth node of the cell membrane is 
calculated using a spring model as follows 

Fbond,f = -^*b-U (3.14) 

where o is the spring constant, l\, is the equilibrium bond length, x\> is the distance spanning the 
endpoints of the GPIba receptor on the platelet surface and the vWF-Al binding site. Table 2 lists 
values of parameters of the DAM used in simulations. 

(c) Lattice Boltzmann method for simulating blood flow 

The LBM uses purely localized fluid particle evolution and relaxation, which in turn facilitates 
parallelization in computer implementation. The LBM decomposes the fluid domain into 
structured lattice nodes and operates on the lattice. The fluid is modelled as a group of fluid 
particles that are only allowed to move between lattice nodes or remain at rest. The composition 
of the lattice nodes depends on the chosen lattice model. In this paper, we used the three- 
dimensional model of a cubic lattice (16 x 64 x 16 |xm with spacing h = 0.2 |xm) with 19 discrete 
velocity directions (model D3Q19, as shown in figure 3). The LBM solves the Boltzmann equation 
describing the dynamics of fluid from a microscopic point of view: in fluid, particles with 
velocities Vj collide with certain probability and exchange momentum. The collisions are assumed 
to be ideal, that is the total momentum and energy is conserved during the collisions. The 
Boltzmann equation describes the probability f{x,v,t) of finding a particle with velocity v at a 
position x and at time t evolving with time 



v-V x / + F.V p /+^ = ^(f), 



(3.15) 




Figure 3. Lattice Boltzmann D3Q19 (three-dimensional and 19 velocities) model. The lattice vectors q represent the velocities 
of the particle moving from the centre grid point to its neighbour grid point. 



Table 2. Values of physical parameters used in simulations. 



I parameters 


definition 


value 


references I 


a 


platelet radius 


1.0 |xm 


[23] 


X 


platelet aspect ratio 


0.25 


[22] 


r 


flow shear rate 


300 and 400 s- 1 


[1] 


p 


blood plasma density 


1.0239 gem" 3 


[24] 


ii 


blood plasma viscosity 


1.2 cP 


[24] 


E 


platelet elastic modulus 


25 kPa 


[11] 


h 


average length of initial spring length 


75 nm 




k s 


global area constraint coefficient 


6000 k^T/l] 


[25] 


k t 


local area constraint coefficient 


moki/i] 


[25] 


kv 


volume constraint coefficient 


6000 k E T/l] 


[25] 


k 


bending modulus 


200 kl 


[27] 


T 


temperature 


300 K 




k° 

"on 


intrinsic cross-linking formation rate 


10-V 1 


[1] 


"off 


unstressed disassociation rate 


3.45 s- 1 


[1] 



where F denotes an external body force, V J/P is the gradient in position and momentum space and 
Q{f) denotes the collision operator which is chosen as a relaxation of / with a characteristic time 
t to the equilibrium distribution/^ (v, p): 

^(/) = "7(/-/ (eq) )- ( 3 - 16 ) 



The equilibrium distribution function depends on the local density p(x, t) and the velocity field 
v(x,t). In the D3Q19 lattice model, 19 values fi(x, t) are stored at each lattice site assigned to a 
lattice vector q. The local density at a lattice point is obtained by summing all//, 



19 



i=l 



(3.17) 



and the streaming velocity is given by 



1 19 

u(x, t) = — — *)q, 



(3.18) 



z=l 



where q = h/ Mis the lattice speed associated with the zth direction and Af is the time step of our 
simulation. 

Using a Chapman Enskog expansion, Guo et al. [37] showed that the following lattice 
Boltzmann equations give a second-order accurate v, the Navier-Stokes velocity in the presence 
of a spatially varying, time-dependent force 



fi(x + a At, t + At) =fi(x, t)-- (/■•(*, t) -ff\p, v)) 



(F-q) (uF T +Fu T ):(cic] + c 2 s I) 



(3.19) 



where u is a streaming velocity defined in equation (3.18), z? = u + FAt/2p, and 



f (eq W)=^P 



C; • » (C, • V) Z 
1+ + 



24 2c? 

with the lattice speed of sound c s = ^1/V3^ h/ At for the D3Q19 lattice and the lattice weights 



(3.20) 



2/36, i = 1 . . . 6, 

^• = 1/36, f = 7 ... 18, 
12/36, z = 19. 



(3.21) 



The pressure p = c^p turns out to be proportional to the density and the dynamic shear viscosity 
is given by 

71 = c 1 s p(x-]\. (3.22) 



To ensure convergence and stability of LBM, we follow the method in [38] to choose our 
parameters. Spacing h = 0.2 |xm was determined by our simulated fluid domain and memory 
size of the GPU card. Time step At was determined from the equation [38]: At = (r — 0.5)/z 2 /(3u), 
where v = /x/p is kinetic viscosity, \x and p are fluid viscosity and density as defined in table 2. 
Generally speaking, a larger value of r leads to a more stable LBM simulation, and r must be 
greater than 0.5. We set r = 1.379 s in our model, such that At — 10 -8 s. 

Periodic boundary conditions in x-z and y-z boundary planes (y = 0, y = 64 |xm, x = 0 and 
x = 16 |xm), are realized by propagating the ft from the computational domain on the one 
boundary to the boundary on the opposite side of the domain. In the x-y boundary planes, we 



used the onsite velocity boundary conditions proposed by Hecht & Harting [39]. For instance, in 
the x-y boundary plane z = 0,fi (i = 1, 2, 3, 4, 6, 7, 8, 10, 11, 12, 14, 16, 18, 19) can be obtained from 
the streaming step, but fj (i = 5, 9, 13, 15, 17) are undetermined. Following the methods of Hecht 
& Harting [39], we obtain 

P = ^—Ifl +fi +h +h +fr +/s +/12 +/19 + 2(/ 6 +/io +/14 +/16 +/l8)L (3-23) 

1 — D z 
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f9=fu + ^(v z + v x )-N z x/ 


(3.25) 


fl3 =flo + ^(v z -v x )+N z x/ 


(3.26) 


fl5 =fl8 + ^(Vz + Vy)-N Z y 


(3.27) 


and fi 7 =/i 6 + ^(u z - uy) + Ny. 


(3.28) 



Here, v x , Vy and v z are boundary velocities in x-, y- and z-directions, and N z the transverse 
momentum corrections on the z-boundary for distributions propagating in x- and y-directions, 
respectively 

K = \lfl +f? +fs - (fi +fu +/12)] - \ P v x (3.29) 

and 

Ny = \[h +^ - (/4 +/ 8 +/12)] - \pVy. (3.30) 

(d) Coupling platelet, dynamic adhesion and flow submodels 

(i) Coupling platelet and dynamic adhesion submodels 

As described in §3c, GPIba receptors are randomly distributed on the cell membrane. In each step 
of the simulation, forming and breaking a GPIba-vWF bond is updated using the DAM. When 
a formed bond is either stretched or compressed, the bond deformation force is computed using 
equation (3.14). This force is exerted on the cell membrane at the place where the GPIba receptor 
of the bond is located. When only the platelet membrane and vessel wall interaction is considered, 
the term H wa \\ of equation (3.1) represents the energy associated with these interactions. In 
particular, 3(H wa n)/3x^ corresponds to the sum of following two forces (i) adhesion forces caused 
by GPIba-vWF bond and (ii) short-range repulsive forces accounting for contact of vessel 
wall. The short-range repulsive force is given by an empirical relationship: F rep = Fo( Te_T£ )/ 
(1 — e~ TS ), where Fq = 500 p Nm, r = 2000 |xm _1 and s is the separation distance between platelet 
membrane and vessel wall [29]. Thus, dH wa \\/dxi term in equation (3.11) is defined as 

^ x — ^bond,/' ~r -'"rep- 

(ii) Coupling cell and flow submodels 

To couple the integrated platelet and stochastic DAM submodels with the blood flow computed 
by the LBM, we use the IBM [16]. In the IBM (figure 4), the Eulerian description is used for the 
fluid dynamics, and the Lagrangian description is used for objects immersed in the fluid. Using 
lowercase letters for Eulerian variables, and uppercase letters for Lagrangian variables, we have 

diX f 

u(x,t)8(x-X)dx (3.31) 



^ :U(X,t)-. 



and 



f(X,t) = 



¥(X,t)8(x-X)dX. (3.32) 

n 



Figure 4. Eulerian fluid grid (black) and Lagrange solid elements (red). A Eulerian description is used for the fluid dynamics, and 
a Lagrangian description is used for objects immersed in the fluid. The communication between these two coordinate systems 
is realized by the immersed boundary method. (Online version in colour.) 



where t is time, u the flow velocity, U the speed of the solid object boundary, x the fluid flow 
coordinate, X the boundary coordinate, / the force density on the fluid node, F the force density 
on the solid elements and 8(r) the Dirac delta function. 

Equations (3.31) and (3.32) are approximated using a regularized discrete delta function 8^. 
The discretized forms of equations (3.31) and (3.32) using 8^ are as follows 

= U k = Yl U ijk 8 h(Xijk ~ Xm)h 3 (3.33) 

i,j,k 

and 

fijk = ¥ ™8 h {xijk ~ Xm)h 3 , (3.34) 

m 

where h is the fluid node spacing, = {ih,jh, kh) the coordinate of the kth Eulerian grid node, 
X m the Lagrange coordinate of the rath elements,/^ the force density on X{p F m the force density 
on X m , Ujjk the velocity of x^, the velocity of X m . The discrete delta function 8^ appearing 
in equations (3.33) and (3.34) is a smoothed approximation to the Dirac delta function 8(r). (The 
detailed derivation procedures in several forms have been presented in the literature [40].) We 
use the following common form 



and 

(p(x) 



Ww) 4*S)*(D*(D (335) 
i( 1+cos (¥))' f -°^ 2 (336) 

0, for \x\ < 2. 



To sum up, first, in each step of the simulation, equation (3.20) is solved. Then, positions of 
nodes of the platelet membrane are updated by equation (3.33). Finally, the MC computations are 
performed to break the already formed bonds and to generate new bonds from unbound GPIba 
and vWF. 

4. Results 



First, the model was verified by comparing simulation results with analytical solutions and 
available model simulation data [24]. Next, we validated the model by comparing the simulation 
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Figure 5. Validation of the platelet dynamics model. The analytical solution for the platelet rotational trajectory (Jeffery 
orbit), trajectory calculated by a completed double-layer-boundary integral equation method (CDL-BIEM) and our simulation 
(LBM-IBM) are shown by solid and dashed colour lines (inset key). (Online version in colour.) 
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Figure 6. The configuration of the simulated platelets at different times flipping over the vessel wall for the wall shear stress 
of 3.0 dyn cm -2 . Image sequence (a) shows projection of the platelet on the x-y plane from dimensionless time point 0 to 
14. Image sequence (b) shows projection of the platelet on they-z plane from dimensionless time point 0-14. The coordinate 
system is defined in figure 2a. (Online version in colour.) 



results with the experimental data [1] on flipping platelets flowing over a vWF-coated surface. 
A calibrated three-dimensional model was used to predict how the stiffness of a platelet 
membrane, the number of receptors on a platelet membrane and the strength of platelet-platelet 
adhesion affect the paused time of the platelet adhering to a vessel wall. 

(a) Validation of fluid-platelet coupling by comparing with the Jeffery orbit 

Mody et al. [10] described theoretical solutions using the Jeffery orbit theory and provided 
predictions obtained using the analytical platelet-flipping model. This analytical solution (shown 
as solid red line in figure 5) did not consider the wall effect and applied only to the cases of 
platelet motion far from the wall (H/a > 20) [24], where H is the centroid height of platelet and 
a is the major radius (as shown in figure 2). Mody et al. [24] modified the completed double- 
layer-boundary integral equation method to include a flat surface boundary that was used to 
compute the effects of the wall on the flow behaviour of a platelet. Platelets located as far as 
2.4-fold platelet radii from the surface 'display modified' Jeffery orbits with periodic rotational 
motion in the direction of flow (green dashed line in figure 5). To verify our model, we simulated 
the flipping of a single platelet located at the distance of 2Aa as well as greater than 20a, from 
the vessel wall. Our simulations revealed that the calculated orbit of rotation (blue dashed line 
in figure 5) agreed well with the Mody's simulation results [24] (green dashed line) within an 
experimental error of 2.65% for platelet located at the distance of 2Aa, and agreed perfectly with 
Mody's simulation results [24] for platelets located at a distance of more than 20a. Figure 6 shows 
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Figure 7. The number of tethering events as a function of the platelet paused time. The solid lines are the fitting lines of 
experimental data for shear stresses of 3.0 dyn cm -2 (shown in blue) and 4.0 dyn cm -2 (shown in red). The corresponding 
slopes of the fits (/r of f values) are —4.83 and —5.18. The dashed lines are the fitting lines of simulation results (shown with 
circles) for shear stresses of 3.0 (shown in blue) and 4.0 dyn cm -2 . The corresponding slopes of the fits (/r 0 ff values) are —3.31 
and —3.58. (Online version in colour.) 
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the series of snapshots from our simulations of platelet-flipping in a shear flow near the vessel 
wall. In our model, the platelet was modelled as an elastic cell with the elastic modulus measured 
by the AFM experiments [11], whereas Mody et ah [10] considered the platelet as a rigid object. By 
comparing our simulation results and results of Mody et ah [10], we conclude that our simulations 
can be successfully implemented to model the motion of individual resting platelets revealing 
high stiffness membrane values. 

(b) Validation of the model of the platelet-substrate adhesion 

To validate the kinetic submodel, we simulated flowing platelets adhering to substrate through 
GPIba-vWF binding and calculated k Q a rates to compare with available experimental data. The 
model parameters used in our simulations (table 2) were obtained in biological experiments 
[1,11,22-25,27]. The adhesive dynamic parameters were measured in in vitro flow chamber 
tests [1]. 

Doggett et al. [1] measured the kinetics that governs platelet interactions with vWF in 
haemodynamic flow. In their experiment, the frequency of tethering for platelets was measured 
by determining the percentage of cells that paused, but did not translocate, on vWF substrates. 
The frequency of tethering for microspheres coated with vWF on antibody-immobilized platelet 
substrates was also measured. A transient tether event was defined as flowing platelets that 
abruptly halted forward motion for a defined period of time and subsequently released, without 
evidence of translocation, to resume a velocity equivalent to that of a non-interacting cell. 
Dissociation rate constants (/c Q ff) were determined by plotting the natural log of the number of 
beads that interacted as a function of pause time after the initiation of tethering (figure 7, the 
slope of the line is — /c G ff). 

To calculate dissociation constants, we performed simulations for various numbers of random 
seeds (1000-1200). The results of the simulations and experimental data are presented in figure 7 



for two different flow shear rates as the natural log of the number of platelet tethering events 
versus the pause time. The values of dissociation rate k 0 # were found to be 3.31 and 3.58 s _1 
for flow shear stress 3.0 and 4.0 dyn cm -2 , respectively. The corresponding experimental values 
obtained in [1] were 4.83 and 5.18s _1 . 

It should be mentioned that our simulations confirmed several experimental observations. It 
was reported in [1] that in the range of flow rates considered, forces acting on the GPIba-vWF 
bond were not sufficient to alter the rate of dissociation k Q ft. Our simulations also demonstrated 
that the k 0 # values altered only in a small range from 3.31 to 3.58 s -1 . Additionally, it was reported 
in [1] that the forces acting on a platelet in shear flow were 14.7 and 19.6 pN for flow shear stress 
3.0 and 4.0 dyncm -2 , respectively, whereas our model yielded very close force values of 12.8 and 
15.6 pN, respectively. 

Our simulations also confirmed that in the range of flow rates (0-4 dyn cm -2 wall shear stress), 
bond association and dissociation kinetics can be successfully described by the Dembo model 
(equations (3.12) and (3.13)). 

(c) Predictive simulations 

The responses of a platelet to interactions with the environment depend, among others, on 
the mechanical forces that platelets experience. Here, we consider effects of platelet membrane 
tension, flow shear stresses and adhesion bond forces on platelet-substrate adhesion dynamics. 

(i) Effect of platelet membrane stiffness 

Simulation results in §4a show that the flow dynamics of the platelet in linear shear flow can 
be studied by modelling platelets as rigid objects. How the stiffness of the platelets affects the 
platelet-substrate interaction remains to be answered. In [41], it was reported that alteration of 
platelet stiffness can modulate platelet aggregation. We hypothesized that softer cells lead to 
prolonged adhesion time and could potentially increase chances of platelets being activated after 
adhesion. Here, we report the simulation results indicating remarkable changes in platelet paused 
time as the platelet membrane stiffness changes. We varied the platelet membrane stiffness 
from 25 to 2.5 kPa, and performed simulations with 30 different random seeds to obtain 30 
different paused times under flow shear stress of 3.0 dyn cm -2 . The paused time was 6.69 ± 0.71 s 
(mean ± SD) for the membrane stiffness of 25 kPa, which is about twice higher than the paused 
time 3.15 ± 0.69 s (mean ± SD) for the membrane stiffness of 2.5 kPa (Mest, p < 0.0008, figure 8). 
The total deviation of all the nodes in the deformed shape in figure 8a is 3.5 |xm compared with 
the reference configuration, and in figure 8b is 0.28 |xm. Thus, these simulation results indicated 
that softer cells have prolonged average paused time. 

(ii) Effect of the number of GPIba receptors expressed on the platelet membrane 

The interaction between platelet glycoprotein (GP) Ib-IX-V complex and vWF is the first step 
of the haemostatic response to vessel injury. As resting platelets interact with vWF, binding 
of vWF to GPIba initiates platelet activation [42]. Meanwhile, in platelet- type von Willebrand 
disease, mutations of GPIb functional receptors can compromise haemostasis by increasing the 
affinity for vWF [43,44]. Some studies demonstrated that abnormalities in the concentrations of 
GPIb membrane proteins are present in patients with myeloproliferative disorders. In particular, 
decreased GPIb concentrations were found in patients with thrombocythaemia and leukaemia 
[45,46]. How the platelet-substrate adhesion dynamics and subsequent platelet activation are 
affected by the number of GPIb is not clear. The objective of our simulations performed in 
this section was to gain insights into this problem. We varied the platelet receptor number 
from 10 688 (normal) to 5344 (insufficient), and performed simulations with 30 different random 
seeds to obtain 30 different paused times under flow shear stress of 3.0 dyncm -2 . The results 
of our simulations revealed that the paused time in the case of decreased receptor number was 
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Figure 8. The simulated deformations of platelet structures during their adhesion to the vessel wall for platelet stiffness of 2.5 
(a) and 25 kPa (b). The effect of the platelet membrane stiffness on the platelet paused time (c). The paused time was 6.69 ± 
0.71 s (mean ± SD) for the membrane stiffness of 25 kPa, which was about twice higher than the paused time of 3.15 ± 0.69 s 
(mean ± SD) for the membrane stiffness of 2.5 kPa. (Online version in colour.) 
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Figure 9. The effect of the number of platelet receptors on the platelet-vessel wall paused time. The platelet paused time for a 
decreased number of GPIb functional receptors was 2.07 ± 0.41 s (mean ± SD), which was significantly lower than the paused 
time of platelets having the normal number of receptors (3.15 ± 0.69 s, mean ± SD). (Online version in colour.) 



2.07 ± 0.41 s (mean ± SD), which was significantly lower than 3.15 ± 0.69 s (mean ± SD) for the 
normal receptor number group (Mest, p < 0.02, figure 9). Our simulations predicted that as the 
number of GPIb on the platelet membrane decreased, the paused time of platelet adhesion to 
vessel wall also decreased. Thus, the results of our model suggest that the number of functional 
GPIb is an important factor determining platelet adhesion and subsequent activation. This has 
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Figure 10. (a) Initial configuration of two platelets used in simulations studying the effect of mutual interaction of platelets on 
platelet-wall adhesion, (b) Platelet-vessel wall paused time as a function of the number of interacting platelets. The platelet 
paused time was 1.61 ± 0.46 s (mean ± SD) for two adhesive platelets, which was significantly lower than the stopping time 
of platelets having the normal number of receptors (3.15 ± 0.69 s, mean ± SD). (Online version in colour.) 



important biological consequences, as controlling the number of functional GPIb receptors can 
provide means for development of novel anti- thrombotic drugs. The mechanism of these drugs 
is based on inhibiting /promoting the function of platelet GPIb receptors to decrease /increase 
adhesion of platelets to vWF to control blood clot growth [47]. 

(iii) Effect of the platelet-platelet adhesion 

To study how platelet-platelet interaction affects platelet adhesion to the blood vessel wall, we 
modelled dynamics of two platelets near the surface of the vessel (figure 10a). In the model, 
the two platelets interacted with each other and one of them adhered to the vessel wall. Our 
simulations revealed that the platelet paused time was 1.61 ± 0.46 s (mean ± SD) in the case 
of two adhesive platelets, which was significantly lower than the pause time of 3.15 ± 0.69 s 
(mean ± SD) calculated for a single platelet interacted with the wall (Mest, p < 0.02, figure 10b). 
These results indicate an important mechanism by which a single platelet adhesion can be affected 
owing to interaction with neighbouring cells. These findings have direct biological consequences 
and help to explain how the increased platelet concentration in blood can affect platelet- 
wall adherence. 

5. Discussion 

This paper described a novel three-dimensional model coupling processes at three biologically 
important spatial scales critical for early blood clot development and uses models to provide 



predictive simulations. First, our model provides a comprehensive representation of mechanical 
properties of a platelet based on the implementation of a hybrid membrane submodel to describe 
mechanical behaviour of the cytoskeleton network and the lipid bilayer of the platelet. In previous 
studies, platelets were modelled as rigid bodies [24,33,48]. However, it has been experimentally 
shown [11] that platelets exhibited both elastic and viscoelastic behaviour and that they undergo 
large deformation in shear flow [12]. 

Experimental studies demonstrated [1,49,50] that flow shear stress could increase both bond 
formation and dissociation rates during platelet adhesion to the vessel wall. Additionally, 
estimates for the forces acting on platelet-substrate bonds were provided in [10]. However, Xu 
et ah [8] did not describe a detailed computational model to simulate the binding dynamics under 
various flow conditions. By combining three-dimensional multi-scale models with microfluidic 
experiments, we provided a methodology to quantify in detail single platelet-flipping in blood 
flow and platelet tethering to the injured vessel wall. It results in a two-way coupled fluid- 
cell interaction submodel combined with a stochastic submodel of formation /breakage of 
individual receptor-ligand bonds. This approach provided a biologically justified description of 
platelet dynamics, which can be used to simulate dynamics of platelets under more complex 
flow conditions. 

By incorporating physiological parameter values characterizing cellular membrane mechanics, 
our method provides explicit representation for the structure of the cytoskeleton and simulation 
of cellular dynamics. Thus, our model allows one to examine how the mobility of cells is affected 
by their membrane structural and mechanical properties and, hence, aids in providing prognostic 
assessment in blood cell disorders outcome. The model developed in this paper can be also 
used for simulating important biomedical problems which involve description of dynamics 
and deformation of cells in fluid flow, including (patho)physiological inflammation involving 
leucocyte and platelet tethering to the vessel wall. Other important applications of the model 
include studying cell aggregate formation in blood, metastasis of tumour cells as well as stem cell 
attachment to the target tissues. 
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Appendix A. GPU implementation 

The CPU used for our simulation is the Intel Xeon CPU L5520 with clock rate 2.27 GHz. Our 
NVIDIA graphics card is the GeForce GTX 480, Clock rate: 1.45 GHz, CUDA driver version: 
5.50, CUDA runtime version: 5.50, CUDA capability version: 2.0. GPUs are separate devices 
with their own processors, and memory devices which do not have direct access to the CPUs 
or CPU memory units. The communication pathway for transferring data between CPU memory 
and GPU memory has a relatively slow bandwidth capability compared with direct access to 
memory devices. Thus, it is necessary to minimize communication as much as possible. The 
typical GPU code is composed of three main parts: (i) initialization, (ii) execution, and (iii) clean- 
up. During initialization, model data are firstly allocated and initialized in the CPU memory. The 
CPU code then initializes connection to GPU device and allocates GPU memory for the model 
simulation data. The model data are copied from CPU memory to GPU memory units. During 
execution, GPU kernel functions are called and occasionally copy data between CPU and GPU 
memory devices. When a simulation is finished, both CPU and GPU memory units are freed and 
connection to GPU device is shut down for the clean-up. 

Figure 1 shows the flow chart of our simulation algorithm. For a single step of the simulation, 
its execution on the GPU starts with the hybrid membrane model. Assuming that a platelet 
consists of N triangle mesh elements and P nodes (each node represents an SCE; figure 2b), The 
GPU kernel function to calculate forces in equation (3.3) has the form 



Table 3. Execution time of 10 000 simulation steps on CPU and GPU for different fluid grid sizes. 
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80 x 300 x 20 
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40 x 150 x 10 
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20 x 75 x 5 


232 


51 



sem_Force_kernel< < < blocks? er Grid, threadsPerBlock >>>(grids->devlmage){ 
//Calculation of forces acting on subcellular elements 



1 

where blocksPerGrid and threadsPerBlock are determined according to the block and thread 
distribution on the GPU card, and P = blocksPerGrid*threadsPerBlock. If this is implemented on a 
single CPU in serial configuration, there will be a total of P iterations for one step of simulation. In 
our GPU implementation, this simulation is performed simultaneously on P GPU threads. Hence, 
it reduces the complexity of execution time from O(P) to O(l) for one step of simulation. Similarly, 
we use the following GPU functions to calculate forces owing to bending, area constraint and 
volume constraint 

sem_Bending_kernel<<< blocksPerGrid, threadsPerBlock >>>(grids->devlmage){ 
//Calculation of forces due to bending 



) 

sem_Area_Volume_kernel<<< blocksPerGrid, threadsPerBlock >>>(grids->devlmage){ 
//Calculation of forces due to area and volume constraint 



) 

where blocksPerGrid and threadsPerBlock are determined according to the block and thread 
distribution on the GPU card, and N = blocksPerGrid*threadsPerBlock. It reduces the complexity 
of execution time from O(N) to O(l) comparing with CPU code. 

The forces acting on solid nodes spread to the neighbour fluid nodes using the IBM. The GPU 
kernel function has the form 

fluid3dJ r orce_distributeJcernel<<<blocksPerGrid3D,threadsPerBlock3D>>>(grids->devImage){ 
//Implementation of IBM 



) 

where both blocksPerGrid3D and threadsPerBlock3D have three-dimensional structure similar to 
a three-dimensional space coordinate. Let blocks per grid 3D = (x^,y^,z^) and threads per block 
3D = (xt,yt,zt)r and the fluid lattice has the size of X x Y x Z, then X = x^ x xt, Y = y^ x yt and 
Z = Zb x Zf. It reduces the complexity of execution time from O(XYZ) to O(l) comparing with 
CPU code. Similar strategy is applied to lattice Boltzmann method implementation. 

There are M receptors in each of the N triangle mesh elements of a platelet. The GPU kernel 
function for dynamic adhesion model has the form: 

sem_platelet_wall_kernel< <<blocksPerGrid, threadsPerBlock> > > (grids->devlmage){ 
//Implementation of DAM for a single receptor 



) 

where N x M = blocks per grid x threads per block. It can speed up the execution time NM times 
compared with CPU implementation. 

In conclusion, the GPU will reduce the whole simulation from O(P) + O(N) + O(XYZ) + 
0(NM) ~ O(XYZ) to O(l) in time complexity. Table 3 shows real execution time of 10000 step 



simulation on both CPU and GPU for three different fluid grid sizes. As table 3 shows, the 
execution time of GPU code is only about 1/100 of the CPU version for the fluid grid size we 
used in this paper. While the CPU execution time increased linearly with fluid grid size, the GPU 
execution time increased much slower than CPU. 

References 

1. Doggett TA, Girdhar G, Lawshe A, Schmidtke DW, Laurenzi IJ, Diamond SL, Diacovo TG. 
2002 Selectin-like kinetics and biomechanics promote rapid platelet adhesion in flow: the 
GPIb(alpha)-vWF tether bond. Biophys. J. 83, 194-205. (doi:10.1016/S0006-3495(02)75161-8) 

2. Hussain MA, Agnihotri A, Siedlecki CA. 2005 AFM imaging of ligand binding to platelet 
integrin alphallbbeta3 receptors reconstituted into planar lipid bilayers. Langmuir 21, 6979- 
6986. (doi:10.1021/la046943h) 

3. Holland NB, Siedlecki CA, Marchant RE. 1999 Intermolecular force mapping of platelet 
surfaces on collagen substrata. /. Biomed. Mater. Res. 45, 167-174. (doi:10.1002/(SICI)1097-4636 
(19990605)45 : 3<167::AID-JBM2>3.0.CO;2-7) 

4. Vinckier A, Semenza G. 1998 Measuring elasticity of biological materials by atomic force 
microscopy. FEBS Lett. 430, 12-16. (doi:10.1016/S0014-5793(98)00592-4) 

5. Yago T et al. 2008 Platelet glycoprotein Ibalpha forms catch bonds with human WT vWF but 
not with type 2B von Willebrand disease vWF. /. Clin. Invest. 118, 3195-3207. 

6. Sun L, Rajamannan NM, Sucosky P. 2013 Defining the role of fluid shear stress in the 
expression of early signaling markers for calcific aortic valve disease. PLoS ONE 8, e84433. 
(doi:10.1371/journal.pone.0084433) 

7. Sun L, Chandra S, Sucosky P. 2012 Ex vivo evidence for the contribution of hemodynamic 
shear stress abnormalities to the early pathogenesis of calcific bicuspid aortic valve disease. 
PLoS ONE 7, e48843. (doi:10.1371/journal.pone.0048843) 

8. Xu Z, Kim O, Kamocka M, Rosen ED, Alber M. 2012 Multiscale models of thrombogenesis. 
Wiley Inter discip. Rev. Syst. Biol. Med. 4, 237-246. (doi:10.1002/wsbm.H60) 

9. Xu Z, Kamocka M, Alber M, Rosen ED. 2011 Computational approaches to studying 
thrombus development. Arterioscler. Thromb. Vase. Biol. 31, 500-505. (doi:10.1161/ATVBAHA. 
110.213397) 

10. Mody NA, Lomakin O, Doggett TA, Diacovo TG, King MR. 2005 Mechanics of 
transient platelet adhesion to von Willebrand factor under flow. Biophys. J. 88, 1432-1443. 
(doi:10.1529/biophysj.l04.047001) 

11. Radmacher M, Fritz M, Kacher CM, Cleveland JP, Hansma PK. 1996 Measuring the 
viscoelastic properties of human platelets with the atomic force microscope. Biophys. J. 70, 
556-567. (doi:10.1016/S0006-3495(96)79602-9) 

12. Lettinga MP, Holmqvist P, Ballesta P, Rogers S, Kleshchanok D, Struth B. 2012 Nonlinear 
behavior of nematic platelet dispersions in shear flow. Phys. Rev. Lett. 109, 246001. (doi:10. 
1103/PhysRevLett.l09.246001) 

13. Sandersius SA, Newman TJ. 2008 Modeling cell rheology with the subcellular element model. 
Phys. Biol. 5, 015002. (doi:10.1088/1478-3975/5/l/015002) 

14. Xu Z, Christley S, Lioi J, Kim O, Harvey C, Sun W, Rosen ED, Alber M. 2012 Multiscale model 
of fibrin accumulation on the blood clot surface and platelet dynamics. Methods Cell Biol. 110, 
367-388. (doi:10.1016/B978-0-12-388403-9.00014-X) 

15. Sweet CR, Chatterjee S, Xu Z, Bisordi K, Rosen ED, Alber M. 2011 Modelling platelet-blood 
flow interaction using the subcellular element Langevin method. /. R. Soc. Interface 8, 1760- 
1771. (doi:10.1098/rsif.2011.0180) 

16. Peskin CS. 1972 Flow patterns around heart valves: numerical method. /. Comput. Phys. 10, 
252. (doi:10.1016/0021-9991(72)90065-4) 

17. Fogelson AL, Guy RD. 2004 Platelet-wall interactions in continuum models of platelet 
thrombosis: formulation and numerical solution. Math. Med. Biol. 21, 293-334. (doi:10.1093/ 
imammb/21.4.293) 

18. Sui Y, Chew YT, Low HT. 2007 A lattice Boltzmann study on the large deformation of red 
blood cells in shear flow. Int. J. Mod. Phys. C 18, 993-1011. (doi:10.1142/S012918310701108X) 

19. Skorczewski T, Erickson LC, Fogelson AL. 2013 Platelet motion near a vessel wall or 
thrombus surface in two-dimensional whole blood simulations. Biophys. J. 104, 1764-1772. 
(doi:10.1016/j.bpj.2013.01.061) 



20. Yago T, Wu J, Wey CD, Klopocki AG, Zhu C, McEver RP. 2004 Catch bonds govern adhesion 
through L-selectin at threshold shear. /. Cell Biol 166, 913-923. (doi:10.1083/jcb.200403144) 

21. Savage B, Saldivar E, Ruggeri ZM. 1996 Initiation of platelet adhesion by arrest onto 
fibrinogen or translocation on von Willebrand factor. Cell 84, 289-297. (doi:10.1016/S0092- 
8674(00)80983-6) 

22. Frojmovic M, Longmire K, van de Ven TG. 1990 Long-range interactions in mammalian 
platelet aggregation. II. The role of platelet pseudopod number and length. Biophys. J. 58, 
309-318. (doi:10.1016/S0006-3495(90)82378-X) 

23. Popel AS, Johnson PC. 2005 Microcirculation and hemorheology. Annu. Rev. Fluid Mech. 37, 
43-69. (doi:10.1146/annurev.fluid.37.042604.133933) 

24. Mody NA, King MR. 2005 Three-dimensional simulations of a platelet-shaped spheroid near 
a wall in shear flow. Phys. Fluids 17, 113302. (doi:10.1063/1.2126937) 

25. Li J, Dao M, Lim CT, Suresh S. 2005 Spectrin-level modeling of the cytoskeleton and 
optical tweezers stretching of the erythrocyte. Biophys. J. 88, 3707-3719. (doi:10.1529/biophysj. 
104.047332) 

26. Wu JS, Aidun CK. 2010 Simulating 3D deformable particle suspensions using lattice 
Boltzmann method with discrete external boundary force. Int. J. Numer. Methods Fluids 62, 
765-783. 

27. Du Q, Liu C, Ryham R, Wang XQ. 2009 Energetic variational approaches in modeling vesicle 
and fluid interactions. Physica D 238, 923-930. (doi:10.1016/j.physd.2009.02.015) 

28. Le DV, White J, Peraire J, Lim KM, Khoo BC. 2009 An implicit immersed boundary 
method for three-dimensional fluid-membrane interactions. /. Comput. Phys. 228, 8427-8445. 
(doi:10.1016/j.jcp.2009.08.018) 

29. King MR, Hammer DA. 2001 Multiparticle adhesive dynamics. Interactions between stably 
rolling cells. Biophys. J. 81, 799-813. (doi:10.1016/S0006-3495(01)75742-6) 

30. Dembo M, Torney DC, Saxman K, Hammer D. 1988 The reaction-limited kinetics of 
membrane-to-surface adhesion and detachment. Proc. R. Soc. Lond. B 234, 55-83. (doi:10.1098/ 
rspb.1988.0038) 

31. Reininger AJ, Heijnen HF, Schumann H, Specht HM, Schramm W, Ruggeri ZM. 2006 
Mechanism of platelet adhesion to von Willebrand factor and microparticle formation under 
high shear stress. Blood 107, 3537-3545. (doi:10.1182/blood-2005-02-0618) 

32. Berndt MC, Shen Y, Dopheide SM, Gardiner EE, Andrews RK. 2001 The vascular biology of 
the glycoprotein Ib-IX-V complex. Thromb. Haemost. 86, 178-188. 

33. Mody NA, King MR. 2008 Platelet adhesive dynamics. Part II: high shear-induced 
transient aggregation via GPIbalpha-vWF-GPIbalpha bridging. Biophys. J. 95, 2556-2574. 
(doi:10.1529/biophysj.l07.128520) 

34. Singh I, Shankaran H, Beauharnois ME, Xiao Z, Alexandridis P, Neelamegham S. 2006 
Solution structure of human von Willebrand factor studied using small angle neutron 
scattering. /. Biol. Chem. 281, 38266-38275. (doi:10.1074/jbc.M607123200) 

35. Fox JE, Aggerbeck LP, Berndt MC. 1988 Structure of the glycoprotein Ib.IX complex from 
platelet membranes. /. Biol. Chem. 263, 4882^890. 

36. Hammer DA, Apte SM. 1992 Simulation of cell rolling and adhesion on surfaces in shear flow: 
general results and analysis of selectin-mediated neutrophil adhesion. Biophys. J. 63, 35-57. 
(doi:10.1016/S0006-3495(92)81577-l) 

37. Guo Z, Zheng C, Shi B. 2002 Discrete lattice effects on the forcing term in the lattice Boltzmann 
method. Phys. Rev. E Stat. Nonliner Soft Matter Phys. 65, 046308. (doi:10.1103/PhysRevE. 
65.046308) 

38. Feng YT, Han K, Owen DRJ. 2007 Coupled lattice Boltzmann method and discrete element 
modelling of particle transport in turbulent fluid flows: Computational issues. Int. J. Numer. 
Methods Eng. 72, 1111-1134. (doi:10.1002/nme.2H4) 

39. Hecht M, Harting J. 2010 Implementation of on-site velocity boundary conditions for D3Q19 
lattice Boltzmann simulations. /. Stat. Mech. (http://arxiv.org/abs/0811.4593) 

40. Peskin CS. 2002 The immersed boundary method. Acta Numer. 11, 479-517. (doi:10.1017/ 
S0962492902000077) 

41. Shamova EV, Gorudko IV, Drozd ES, Chizhik SA, Martinovich GG, Cherenkevich 
SN, Timoshenko AV. 2011 Redox regulation of morphology, cell stiffness, and lectin- 
induced aggregation of human platelets. Eur. Biophys. J. 40, 195-208. (doi:10.1007/ 
S00249-010-0639-2) 



42. Kroll MH, Harris TS, Moake JL, Handin RI, Schafer AI. 1991 von Willebrand factor 
binding to platelet Gplb initiates signals for platelet activation. /. Clin. Invest. 88, 1568-1573. 
(doi:10.1172/JCI115468) 

43. Kumar RA, Dong JF, Thaggard JA, Cruz MA, Lopez JA, Mclntire LV. 2003 Kinetics of 
GPIbalpha-vWF-Al tether bond under flow: effect of GPIbalpha mutations on the association 
and dissociation rates. Biophys. J. 85, 4099-4109. (doi:10.1016/S0006-3495(03)74822-X) 

44. Andrews RK, Lopez JA, Berndt MC. 1997 Molecular mechanisms of platelet adhesion and 
activation. Int. J. Biochem. Cell Biol. 29, 91-105. (doi:10.1016/S1357-2725(96)00122-7) 

45. Mazzucato M, De Marco L, De Angelis V, De Roia D, Bizzaro N, Casonato A. 1989 Platelet 
membrane abnormalities in myeloproliferative disorders: decrease in glycoproteins lb and 
lib /Ilia complex is associated with deficient receptor function. Br. J. Haematol. 73, 369-374. 
(doi:10.1111/j.l365-2141.1989.tb07755.x) 

46. Jensen MK, de Nully Brown P, Lund BV, Nielsen OJ, Hasselbalch HC. 2000 Increased 
platelet activation and abnormal membrane glycoprotein content and redistribution in 
myeloproliferative disorders. Br. J. Haematol. 110, 116-124. (doi:10.1046/j. 1365-2141.2000. 
02030.x) 

47. Ji X, Hou M. 2011 Novel agents for anti-platelet therapy. /. Hematol. Oncol. 4, 44. (doi:10.1186/ 
1756-8722-4-44) 

48. Mody NA, King MR. 2008 Platelet adhesive dynamics. Part I: characterization of platelet 
hydrodynamic collisions and wall effects. Biophys. J. 95, 2539-2555. (doi:10.1529/biophysj. 
107.127670) 

49. Alevriadou BR, Moake JL, Turner NA, Ruggeri ZM, Folie BJ, Phillips MD, Schreiber AB, 
Hrinda ME, Mclntire LV. 1993 Real-time analysis of shear-dependent thrombus formation and 
its blockade by inhibitors of von Willebrand factor binding to platelets. Blood 81, 1263-1276. 

50. Kroll MH, Heliums JD, Mclntire LV, Schafer AI, Moake JL. 1996 Platelets and shear stress. 
Blood 88, 1525-1541. 



